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Abstract 

We investigate particle production from coherent oscillation by using the method 
based on the Bogolyubov transformation. Especially, we study the case when the 
amplitude of the oscillation and also the coupling constants with the oscillating field 
are small in order to avoid the non-perturbative corrections from the broad parametric 
resonance. We derive the expressions for the distribution functions and the number 
densities of produced particles at the leading order of coupling constant. It is, however, 
found that these results fail to describe the exact particle production eventually due to 
the non-perturbative effects even if the coupling constants are small. We then introduce 
a simple method to handle with such corrections, i.e., the time averaging method. It is 
shown that this method successfully provides the evolution of the occupation numbers 
of the growing mode. Further, we point out that the approximate results by this 
method satisfy the exact scaling properties coming from the periodicity of the coherent 
oscillation. 
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1. Introduction 



Coherent oscillation of scalar field plays an important role to describe various phenomena 
in particle physics and particle cosmology. One of the most important examples is the so- 
called slow-roll model of the inflationary universe .®'® A scalar field, called as inflaton, 
is initially displaced from the potential minimum and its vacuum energy leads to the de- 
Sitter expansion of the universe. This class of models elegantly solves the flatness and 
horizon problems of the standard Big Bang cosmology. Moreover, it can give an origin of 
the density fluctuations which is strongly supported from the recent measurements of the 
cosmic microwave background radiation.® After the inflation ends, the inflaton starts to 
cause coherent oscillation around its potential minimum. 

The energy of coherent oscillation is diluted due to the expansion of the universe as well as 
the energy transfer to particles though interaction of the oscillating field. Produced particles 
are then thermalized and the hot universe can be realized. The whole of these processes is 
called as the reheating. In particular, the reheating of the slow-roll inflation gives an initial 
condition of the standard Big Bang cosmology. Therefore, the reheating process is crucial 
for understanding the very early universe. 

In this paper, we focus on the first stage of the reheating, i.e., particle production from 
coherent oscillation. This process has been widely discussed based on that coherent oscilla- 
tion is considered as non-relativistic scalar particles,® and particles are produced through 
their decay and/or scattering processes.® In this case the number of parent non-relativistic 
particles is given by the energy of coherent oscillation divided by the mass of the oscillating 
field, and the number of produced particles is determined by it. On the other hand, it has 
been pointed out that, when the coupling with produced particles and also the amplitude 
of the oscillation become large, the non-perturbative effect becomes significant in the early 
stage of particle production.®"® This process is called as the preheating.® Especially, the 
explosive production of bosonic degrees of freedom can happen due to the broad parametric 
resonance effect. On the other hand, the fermion production at the preheating has been also 
investigated. 113 "^ 

The purpose of this paper is to investigate the production of scalars and fermions from 
coherent oscillation, especially when the coupling constants of oscillating field are very small 
to avoid the effect of the broad parametric resonance. For this purpose, we apply the method 
based on the Bogolyubov transformations.^'^® In this case, the equation of motion for the 
mode functions of produced particles in the presence of coherent oscillation is solved and the 
growth of the mode functions are then interpreted as the production of particles. First, we 
will present the analytical formulae for the distribution functions and the number densities 
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of produced particles by using the perturbative expansion of the coupling constants. We will 
also discuss the conditions under which the perturbative results are justified. Indeed, it will 
be shown that the leading-order results collapse in the end. 

This is a signal that the non-perturbative effect becomes important even if the cou- 
pling constant is sufficiently small. Such a correction is crucial for describing the statistical 
properties of produced particles, namely the effects of the Bose condensation for the scalar 
production and the Pauli blocking for the fermion production. In order to handle the an- 
noying non-perturbative effects we will present the time averaging method, which is familiar 
in the nonlinear dynamical system.^ It will be demonstrated that this method is powerful 
to extract the characteristic evolution of the occupation number for the growing mode, i.e., 
the exponential growth for scalar production or the oscillation between to 1 for fermion 
production. Furthermore, we will show that the results by the time averaging method obey 
the exact scaling property, which is obtained from the periodicity of the the equation of 
motionP This gives a justification for the use of the time averaging method. Throughout 
the present analysis we neglect the expansion of the universe for simplicity. 

The rest of this paper is organized as follows. In Sec. H] we explain the model in this 
analysis. We perform in Sec. E]the perturbative estimation of the yields when the amplitude 
of the coherent oscillation is sufficiently small. The importance of non-perturbative effects 
in particle production is addressed in Sec. HI We present the time averaging method to 
deal with such effects, and try to figure out the statistical properties of produced particles. 
Finally, the last section is devoted to conclusion. We also add Appendix [A] to explain the 
perturbative estimation of the number density. 

§2. Framework 

To begin with, let us explain the framework of this analysis. We shall study the pro- 
duction of real scalar field x an d Dirac fermion ip from the coherently oscillating by using 
Lagrangian 

£ = \ (<V>)(<9» - v(<f>) + l - (d.xWx) -\gl<P 2 x 2 + i^l^ -g F MiJ>, (i) 

where g$ and g F are coupling constants. For definiteness, we take here the potential for the 
real scalar field <fi as 

V{4>)= l -ml{^-(4>))\ (2) 

where m<$, and (0) are mass and vacuum expectation value (vev) of 0, respectively, and they 
are taken to be real and positive. At the potential minimum x an d ip receive masses as 
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m x = gs(4>) and = gF (<!>)■ From now on the field is assumed to oscillate coherently 
around (0) with an amplitude $ 



<f>(t) = (</>)+$ cos(rV) , (3) 

and it is treated as a classical background field. Notice that we neglect the expansion of the 
universe throughout the present analysis. 

Particle production from the oscillation is usually discussed as follows: The coherent 
oscillation is considered as non-relativistic particles.® In this case, the energy density of the 
oscillation is p^ = 4> 2 /2 + V(<p) = m 2 ^ <P 2 /2 (here and hereafter the dot denotes a derivative 
with respect to time), and the number density of is estimated as = p^jm^ = m ( f > (p 2 /2. 
Decays of 0, — > \ + x an d — > ip + ip, are important processes to produce x an d ip- The 
partial rates of these processes are found from Eq. ([T]) as 

1 ^x+x — ~T> > l 4 J 

87T m-A 



9% Pi 

F^+Q = , (5) 



57T 



where f3 x ^ = J 1 — 4m^ ^/m|. The number densities of x an d ip from the decays of are 
then estimated as 

g 2 ^ f3 *p 2 tti 2 

"xW = 2 Wx B ^= S X gn (6) 

n^t) = 2 t = 1- 1 . (7) 

It is seen that proportional to i by neglecting the decrease of rz^, and that they are 

induced at the order of coupling squared. 

Further, the scattering processes of 0's are another sources of particle production. The 
number density of x due to the process + 0— )• x + X is estimated as n x (t) = (cv^) n^t, 
where is the relative velocity of 0's and (crv^) is the invariant scattering rate which is given 
by (avj,) = <?| (3' x /(32 tt m|) with f5 x = yjl — m^/m^ in the non-relativistic limit — > 0. 
We then find that 

gt (5i$ A 

n ( t ) = ^-^ — t, (8) 
xW 128 7T ' K J 

which is again proportional to t, while it is induced at the fourth order of coupling constant. 
Thus, the production via scattering can be neglected as long as the oscillation amplitude is 
sufficiently small, say <C (0). It should be noted that the scattering rate of + 0— > ip + ip 
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vanishes in the non-relativistic limit, and hence the production of ip via scattering is less 
significant. 

It has been discussed in the literature that particle production from the coherent os- 
cillation is more involved than the above naive treatment. In the following, we study the 
production of x an d ip by using the method based on the Bogolyubov transformation.^'^ 
Especially, we concentrate on the case in which the coupling constants g$ and gp are very 
small, say gs,F rn^/'P , in order to avoid the non-perturbative effect due to the broad 
parametric resonance. We perform both analytical and numerical estimations of the yields, 
and find the validity of the naive argument of particle production. 

§3. Perturbative Estimate of Yields 

We are now at the position to derive the analytical expressions for the yields of \ and 
if) at the leading order of the coupling constant g$ or gp. Hereafter, we identify the masses 
m x and rrif as parameters being independent on gs and gp, although they are O(gs) and 
O(gp) quantities. Moreover, we assume that the amplitude of coherent oscillation is small 
as <P ((f)) in order to avoid the production from the scattering processes. *3 

3.1. Production of scalar 

Let us first consider the production of the scalar x- m the presence of coherent oscillation 
in Eq. the equation of motion for x is given by 

p-M*(t)] X (t,x) = 0, (9) 

where M x (t) = gs4>if) = m x + 9s$ cosijn^t) . To solve Eq. (JSJ) we expand \ as 

X(t,x)= j^-^ X k{t)a{k) + xl{t)a\-k) e **, (10) 



(2tt) 3 / 2 

where a{k) and a)(k) are annihilation and creation operators, respectively, and they satisfy 
the commutation relation [a(ki) , a^(k 2 )] = <5 3 (fci — k 2 ). The mode function Xki obeys the 
following equation of motion 

x k (t)+4(t) Xk (t) = o. (ii) 

Here the time-dependent frequency is given by 

ul(t) = k 2 + M 2 (t) = u 2 + 2g s <Pm x cos(m^t) + g 2 s <P 2 cos 2 (m t) , (12) 



*) The case of the large amplitude 3> ((/)) will be discussed in elsewhere.^ 
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where u x = k 2 + m x and k = \k\. In the method based on the Bogolyubov transformation, 
the growth of the mode function corresponds to creations of x'sJ-^ 1 Indeed, the phase-space 
distribution function of produced x's is given by (see, e.g., Ref. [7J)) 

/*(*, k) = \ \x k (t)\ 2 + u 2 k (t) \ Xk (t)\ 2 } - \ . (13) 

The number density of x is then estimated as 

d 3 k 

As the initial conditions we take a plain wave solution such that 



"x(*) = / 7S3s/x(*.*)- ( 14 ) 



^^ = ^=7^' X*(0) = -<w fc (0) Xfc (0). (15) 

In this case, / x (0, /c) = for all the momentum and the initial abundance of \ is zero. 

Now we estimate the distribution function and the number density at the leading order 
of gs- For this purpose, we rewrite \k in the form (see, e.g., Ref. [7])) 

V^(t) v 7 ^! 



;i7) 



It is then found from Eq. f lTTj) that a k and satisfy the equations 

zu k [t) 

fo(t) = P^e-»ti*'*Ma k (t), (18) 
zu k [t) 

where a k (0) = 1 and /3 k (0) = 0. The coefficient functions obey the normalization condition 
!«£:(£) 1 2 — \ fik(t)\ 2 = 1- In this case the distribution function f x is written in terms of (3 k as 

f x (t,k) = \(3 k (t)\ 2 - (19) 
It should be noted that the factor Co k /{2uj k ) in Eqs. (JTTJ) and ([I8~]) is O(gs): 

+ 0(g 2 s ) . (20) 



The initial conditions of and then show that the leading contribution to (3 k is O(gs), 
and the time-dependence of a k appears at 0(g 2 s ). Therefore, f3 k at the leading order is 
obtained as 



p k ( t ) = i ^™* m * f dtx [ e ^- 2 ^ - e -K>+2-x)*i] + 0{g 2 s ) . (21) 



x 
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Fig. 1. Evolution of the occupation number 
f x for the growing mode with k = k*. The 
red solid line shows the numerical result 
while the blue dashed line shows the ana- 
lytical result ([23]) . Here we take gs = 10~ 8 
and f3 x = 0.5. 




Fig. 2. Distribution function f x in momen- 
tum space {k/k*) at z osc = 10 (the red 
solid line) and at z osc = 50 (the blue 
dashed line). Here we take gs = 10 -8 and 



0.5. 



This clearly shows that (3k, or equivalently f x , causes oscillation for all the modes, except 



for the mode with u x = i.e., with the momentum k = fc* where 



(22) 



if m x < 771^/2. *3 In this case the imaginary part of (3k grows linearly with time due to the 
cancellation of the phases between the <fi oscillation and the frequency of a pair of x- It is also 
important to note that u x = for the growing mode suggests the energy conservation 

in the process that <fi at rest decays into a pair of x m the true vacuum. For the growing 
mode the leading contribution to the occupation number is then estimated from Eq. (|T9|) as 



f x (t, k*) 



(23) 



where we have neglected the oscillation terms. 

In order to confirm the obtained results we numerically solve the equation of motion (fTTj) 
and estimate the yield of x- As representative values we will take = 1.5 x 10 13 GeV and 
<P = 3.4 x 10 19 GeV from now on. In Fig. [1] we show the evolution of the occupation number 
for the growing mode with k = k* in terms of the number of 4> oscillation z osc = m^t/(27r) by 



taking g s = 10 8 and (3 X = 0.5 



i.e., m x ~ 



0.43 m^). It is seen that the analytical estimation 



*) For m x ^> rn^, there is no growing mode. The study for such a case will be done in elsewhere.^ 
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in Eq. (123]) successfully abstracts the characteristic behaviour of the occupation number 
f x {t,K) oc t 2 . 

The exact expression for the distribution function at 0(g 2 s ) is obtained from Eqs. ffl~9] 
and (12"T1) . which is given by 



g 2 s <P 2 m 2 , m 2 
16 w 4 



fx(t, k) = —i x I(t, u x , m*) , (24) 



where 

I(t, u x , m ) = - 2 _ 2 J 3mJ + 4c^ + (mj - 4w^) cos(2m^) 

(m^ 4U x j I 



-2m„ 



(rricf, + 2u x ) cos((m^ — 2u x )tj + (m^ — 2w x ) cos((m0 + 2uj x )£) >. (25) 



Fig. [2] shows the numerical results of the distribution functions at z osc = 10 and 50. We have 
confirmed that Eq. (T24"]) agrees with the numerical results for the parameter choice in the 
figure. It is seen that f x has a peak at k ~ fc*. For the case of z osc = 10 the peak is located 
at the momentum slightly smaller than k* because of the effect of the oscillation terms, and 
such an effect become negligible for larger 2; osc . It is interesting to note that the modes with 
k > are produced, which are kinematically forbidden in the process where <j) at rest decays 
into a pair of %■ However, their occupation numbers are highly suppressed and it scales as 
A;~ 8 or k~ 6 when 2z osc is integer or not, respectively. On the other hand, for the modes with 
k k* the occupation number is independent on k and they oscillate around a constant 
value. Furthermore, we can also see from Fig. [2] that the typical width of the peak in f x is 
inversely proportional to time. 

The number density is then found at the leading order 

g 2 $ 2 m 2 f3 x 

n x(t)- — *• ( 26 ) 

Here we have listed only the terms proportional to t and neglected the oscillation terms. 
The derivation of Eq. (126]) is explained in Appendix [Al In Fig. [3] we compare Eq. (126]) with 
the exact numerical result. First, we observe that the number density grows as t 4 initially, 
and there is a descrepancy between the leading 0(g$) and numerical results. After a few 
coherent oscillations, however, the number density is approaching to the 0{g 2 s ) result ( I26p . 
and linearly grows in time with a small correction of oscillation. Moreover, it should be 
noted that the perturbative result (126]) coincides with the naive result given in Eq. (]6]), and 
hence n x at the leading order can be estimated by the decays of non-relativistic into pairs 
of x after a few oscillation. The same conclusion has been obtained in study of the narrow 
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Fig. 3. Evolution of the number density n x in terms of z osc . The red solid line shows the numerical 
result while the blue dashed line shows the analytical result (|26p . Here we take gs = 10 -8 and 
P x = 0.5. 

parametric resonance by using the method of the density matrix.^ Finally, Fig. [3] also 
shows that the perturbative result breaks down for z osc > 100. Therefore, we have clarified 
that the perturbative estimation (as well as the naive estimation in Sec. |2]) of the yield can 
be applied only in the limited time interval. 

3.2. Production of fermion 

Next, we turn to consider the production of the Dirac fermion ip. The equation of motion 
for ijj is 

i 7 '*9 M -M^(*)]v(t,f) = 0, (27) 

where M^(t) = gp(j){t) = + gp^> cos(m^t). We decompose ip as 

d 3 k 



4>(t,x) 



(2tt) 3 /2 



Wit, k) b h (k) + v h (t, k) d\{-k) 



ik-x 



(28) 



where the summation is taken over helicity h — ±, and Vh{t,k) = u c h {t, —k). bh and dh are 
annihilation operators of particle and anti-particle, respectively, and they satisfy the anti- 
commutation relations {^(fci), b\ 2 {k,2)} = {^(^i), ^(^2)} = Sh^^i^i — k 2 ). We write 
the wave function Uh as 



u h (t, k) 



Ph(t, k) 

Q h (t, k) 



® <Ph{k) 



(29) 



where the helicity eigenfunction satisfies kiph(k) = hkifh(k). Note that the mode functions 
Ph and Qh obey the normalization condition \P h (t,k)\ 2 + \Qh(t,k)\ 2 = 1. It is then found 
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from Eq. (1271) that the equation of motion Ph is 



P h (t,k)+u 2 k (t)P h (t,k) = 0. 



where u k is given by 



^ifc(*) = u k(t) + i ^(t) = k 2 + M 2 (t) + i M^t) 



(30) 



(31) 



As in the case of the scalar production, we assume a plain wave solution initially and impose 
the conditions 



P h (0,k) 



lcu k (0) + M40) 
2oo k (0) 



P h (0,k) = -iu k {0)P h {0,k). 



(32) 



It is then seen that the equations of motion and the initial conditions for the h — ± states 
are the same, and so P+(t, k) = P_(i, k). On the other hand, Qh is obtained from P h as 



Qh(t,k) = -i£ PiAt.k) -f iM,-(t) P h [t.k) 



(33) 



This means that Q + (t,k) = — Q-(t, k). The distribution function of ip can be written in 
terms of the mode functions as (see, e.g., Ref. ITS])) 



2u k (t) 



(34) 



where f?^ is defined by 



O h (t,k) = -hk P h (t,k)Q* h (t,k) + Q h (t,k)PZ(t,k) +M i ,(t) \P h (t,k)\-\Q h (t,k) 



-21m 



P h *{t,k)P h (t,k) -M^t) 



(35) 



This clearly shows that the distribution functions for two helicity states are exactly the same. 
Finally, the number density of if) is given by 

d 3 k 



nj,(t) = 4 



(2tt) : 



U(t,k) . 



(36) 



Here a factor of four counts the number of internal degrees of freedom of ijj. 

We then turn to estimate the leading contribution to the yield of ip when gp is very 
small. Since P + — P_, the index h will be implicit from now on. In order to solve Eq. (130]) 
we express P as 



(37) 
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In this case the coefficients A k and B k obey the coupled equations 



Bk{t) = 2^iT) exp 

and their initial values are found from Eq. ( 132|) as 



+2i / dtiUk(ti 



-2i / dtxUkih) 



A fc (0) = + + (wj + 2g F <Pm i> + g 2 F <P 2 ) 



B k (t) , 
A k {t), 

1/2 



(38) 
(39) 



Sfc(O) = 



(40) 



where uj\ = k 2 + m^, is the frequency of if) in the true vacuum <& = 0. From now on let us 
find solutions of A k and B k in power series of coupling g F as we did in the x production. 
The leading term of B k is found to be 0(g F ), which can be obtained as 



B k (t) 



2u k (ti) 



g F $rn<p ,(o) 



A^lf(t,k) + 0(g 2 F ), (41) 



where A^ = y/oj^ + denotes the O{g F ) term of A k which is independent on time, and 
we have introduced 



lf(t,k) 



(42) 



Then, we can again see that B k s oscillate with time for all the modes expect for the mode 
with = m ( f ) /2, i.e., with the momentum 



2 1 



(43) 



if m^p < rrifp/2. Notice that it corresponds to the momentum of ip in the decay process 
— > if) + if). For this growing mode we find that 



QF^A^ 

2 77^ 



(44) 



and increases linearly in time (with corrections of oscillations), which is consequence of the 
cancellation of phases between the (f> oscillation and the frequency of a pair of if). Moreover, 
A kt increases as t 2 at 0(g F ), which can be seen as follows: We find from Eq. (138]) that 



A k (t) = A k (0) + / dt 



u k (t 



2u k (t 



A k (0)-i 9 -4^lf(t,k) + 



■>ooz 



(45) 
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Fig. 4. Evolution of the occupation number 
L, for the growing mode with k = k*. The 
red solid line shows the numerical result 
while the blue dashed line shows the ana- 
lytical result (|4"8|) . Here we take qf = 10~ 8 
and /3 X = 0.5. 




Fig. 5. Distribution function /,/, in momen- 
tum space (k/k*) at z osc = 10 (the red 
solid line) and at z osc = 50 (the blue 
dashed line). Here we take gF = 10 -8 and 



0.5. 



where 

Jj(f, k) = I dt x \(m^ - 2m<p) e *("V+ 2 ^)*i + ( m<f) + 2m^) e -*(^- 2 <^)*il B^) . (46) 
Jo >- 

By using Eq. (j4"4"|) is given by 

A k M = A kt (0) - UF 7 k * t 2 + -... (47) 

o 

Therefore, the leading O(gp) contribution to the occupation number of the mode k = k* is 
found from Eqs. (1541) and (1551) 

MUW^Y* 2 . (48) 



by neglecting the oscillation terms. 

We then compare the above result with the numerical solution of Eq. (130]) which includes 
the higher order terms of gp and the oscillation terms. We show in Fig. H] the evolution 
of the occupation number for the growing mode by taking gp = 10 -8 and (3^ = 0.5 (i.e., 
m^p ~ 0.43 m^). We can see that the perturbative result ( 148|) gives a good approximation 
for fy, and it increases at t 2 with corrections from the oscillation terms. 

The exact expression for the distribution function at the leading gp order is found from 
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Fig. 6. Evolution of the number density nw, in terms of z osc . The red solid line shows the numerical 
result while the blue dashed line shows the analytical result (|50p . Here we take gp = 10 -8 and 
A/> = 0.5. 



Eqs. (HU) and (35J) as 



U{t, k) 



g 2 F <t> 2 rr& fc 2 
16 w} 



(49) 



where I(t, Ufy, m^) is given by Eq. (j251) . Fig. [S] shows the numerical estimation of the dis- 
tribution function /^,. Notice that we have confirmed that the analytical result in Eq. (|4"9"|) 
agrees with the numerical one, as in the case of the scalar production. It is seen that has 
a peak at k ~ and its height scales as t 2 as expected. Further, the figure shows that the 
typical width of the peak becomes narrow as Ak(t) oc 1/t. It is also interesting to compare 
Eq. (T24"|) with Eq. f l4~9~]) . The differences between and f x are in the prefactor {k <B- m x ) and 
in the arguments of the function /. Because of these differences, for the modes k <C k» 
does dependent on k in contrast to scalar production and the suppression of for k ^> k* 
is relaxed. Moreover, the former difference provides the additional (3 2 factor for the number 
density of produced ip. As shown in Appendix |A] the leading contribution to the number 
density of ip is given by 



g 2 F ffi<P 2 m 2 



S7T 



(50) 



by neglecting the oscillation terms. It is then found that the number density becomes 
independent on for <C (as long as <C (0)). Notice that it coincides with the 
naive result in Eq. (J2J). 

In Fig. O the evolution of the number density is also shown. We can see that the 
numerical result is approaching to the analytical estimate f )50|) after a few oscillations of 0. 
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Therefore, the leading 0(g"p) contribution to the if) yield is also described by the decays of 
non-relativistic particles into pairs of if) and if). However, it is also found from Fig. [6] that 
the 0(g 2 F ) result in Eq. ( 150|) breaks down in sufficiently later times, say z osc > 100, as in the 
scalar production. This issue will be discussed in the next section. 

§4. Non-perturbative Corrections 

We have so far estimated the leading contributions to the yields of \ and if) from the 
coherent oscillation. When the amplitude of the oscillation is small enough, the yields are 
induced at the order of 0(g 2 SF ). In this case there exists the growing mode with k = k* 
if m x ^ < and its occupation number grows at the rate t 2 after a small number of 

oscillations. The number density n x ^ becomes proportional to t, which is consistent with 
the naive estimation in Eqs. fl6]) and (JZJ) based on the decays of non-relativistic particles. It 
is, however, found that the perturbative result eventually fails to describe the evolution of 
the yield at later times. 

As for the x production, Fig. [7] shows how the occupation number for the growing mode 
evolves at later epoch. We can see that f x (t, A;*) starts to grow exponentially, which is differ- 
ent from the O(g^) result given in As already pointed out in Ref. [7j), this discrepancy 
comes from the non-perturbative effect due to the narrow resonance at the preheating stage. 
Such an explosive production is reflected the statistical property of x, i.e., the effect of the 
Bose condensation. Then, the number density of x a ls° grows exponentially as also shown in 
Fig. [31 On the other hand, with regard to the if) production, the later evolution of f^(t, fe*) 
is shown in Fig. [HJ In this case the 0(g 2 F ) result in (j4"5j) becomes inconsistent in the end and 
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it oscillates between zero and unity. This oscillation behaviour had already been observed in 
Ref. IT2]) . which is a consequence of the Pauli-blocking effect, i.e., is forbidden to exceed 
unity Accordingly, the number density stops to grow at some point as also shown in Fig. [6j 
The importance of these non-perturbative corrections have already been addressed in 
various cases of the preheating process. It should be stressed that such corrections become 
significant eventually even if the coupling constant gs,F is extremely small. In general, it is a 
difficult task to derive the analytical expression for the yield including the non-perturtative 
effect. The previous works^®'^ had investigated by using the knowledge of the Mathieu 
function,^ since the solution of ( fTTj) is approximately given by this function. Here we utilize 
another method, the time averaging method, which is familiar in the study of the non- 
linear differential equations)-® From now on we will demonstrate that the evolution of the 
occupation number for the growing mode can be successfully described by this method. Note 
that the time averaging method has already been used to describe the resonance structure 
of Mathieu function^ where the authors have focused only on the scalar production, and 
have obtained the Mathieu characteristic exponent of the first instability band and Eq. (I57p 
shown in below. 

We will demonstrate below a simpler application of the time averaging method together 
with the method of variation of parameters. This allows us to clearly derive the analytical 
forms of the mode function and the occupation number for the growing mode. Further- 
more, our approach is applicable not only to the scalar production but also to the fermion 
production as shown in below. 

4.1. Growing mode of scalar 

Let us recall the equation of motion for Xk (TTTj) for the growing mode k = k* 



where we have introduced r = and q x = 2g s < ^/m^. We shall solve this equation 

by using the methods of variation of parameters and time averaging. For this purpose, we 
introduce u\ and U2 by 



dr 2 



Xk, 



+ [l + 2g x y/l-0* 



cos(2r) + q 2 cos 2 (2r) Xk. = , 



(51) 




(52) 



with the condition 




(53) 
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where y\ = cost and y 2 = sinr are the solutions for Eq. ( l5Tj) with q x = 0. We then obtain 
the equations for u\ and w 2 as 




2q x Jl-(5l cos(2r) + q 2 cos 2 (2r) 



smr cos t 

— COS 2 T 



sin 2 t 



— sm r cos r 




(54) 



Now we are interested in the characteristic behaviour of Xk, due to the non-perturbative 
effect and, as shown below, its typical time scale is given by r ~ l/g x which is much longer 
than the period of the <fr oscillation in the weak coupling limit (say, q x <§C 1). In this situation 
we can apply the time averaging method which extract the underlying behaviour over a long 
time scale by integrating out the effects of the rapid oscillations. The u\ and u 2 averaged 
over the oscillation period are denoted by U\ and u 2 , and they satisfy 



(55) 



Here and hereafter, we neglect the contributions of higher order of q x . By using the initial 
conditions ( TT5]) at O(gg), xk after the time averaging is obtained as 





XkA T ) = 2/i( r )wi(r) +y2(r)u 2 {T) 
1 



cosh 



0x< 



ft 



sinh 



'1 



$1 



T 



2 J V 2 

Finally, Eq.f JT3|) gives the occupation number for the growing mode with k = fc* 

'9s® r, ^i(9s$m x 



(56) 



as 



f x \t, 



sinh 



I -Pit 



sinh 



t 



(57) 



2 v «• y \ 

We have checked that the obtained result can describe the exact numerical one in 
Fig. [7] apart from tiny corrections of oscillation. Interestingly, this correctly reproduces 
the initial behaviour in Eq. (T2"5|) for <C t <^ t cr = m^/^gs^m^). When g s = 10~ 8 



and m x ~ 



0.43771^ shown in Fig. [7J the critical time is estimated to be z c 



160 



i.e.. 



^cr — 



7 x 10 GeV~ ). Note that the exact numerical estimation of f x shows f x oc t 
for the beginnings of production (within one oscillation), which can not be explained by 
this result. On the other hand, for t ^> t cr , the occupation number grows exponentially as 
exp(2gs'Pm x t/m l j > ). This exponent is consistent with the result from the narrow paramet- 
ric resonance)-^ Notice again that such non-perturbative correction becomes significant for 
t > ^ cr even if the coupling g s is extremely small. Accordingly, the number density starts to 
grow exponentially at t ~ t cr . 
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It is important to note that the mode function for the growing mode in Eq. (1561) satisfies 
the exact scaling property^ 1 in its time evolution, which is a consequence of the periodicity 
of the cf) oscillation. To see this point, let us first recall the exact equation of motion for Xk 
in Eq. ( fTTj) and denote by Xk\t) an d xf\t) its two linearly independent solutions with the 
initial conditions xi^O) = xf\fy = 1 and x^O) = xi (0) = 0. In this case, Xfc can be 
written without loss of generality as 

Xk(t) = Xk(0) [xP® - luj^xfit)] . (58) 
Due to the periodicity of ujk(t), the independent solutions satisfy the exact scaling property^ 1 




(59) 



Xi\t + T) \( x[ l \T) xi 1} (T) 
xt\t + T) ) \ xt\T) X?CO 
where T is an oscillation period of 4>(t) (i.e., Uk{t + T) = u>k{t)). This shows that we can 
extrapolate the mode function at t — nT by using the solution at t — T recursively. From 
this exact property, the occupation number of x a t the time t = nT is written by 



UnT, k) 



( sinh (nD) 



4u 2 k {T) V sinh (D) 
sinh (nD) 
sinh ((n — 1) D) 



xi\T)+ul(T)xt\T) 
f x ((n-l)T,k), 



(60) 



where D = cosh _1 (x[. 1 ^(T)). Therefore, the occupation number at t — nT can be obtained 
by using the mode function at t = T. 

Now, we examine whether the occupation number for the growing mode by the time 
averaging method ([55]) satisfy these scaling properties or not. It should be noted that our 
result can be written as Eq.(l58"|) with 



*£?(*) = cosh 



g s @m 



x. 



t cos 



sinh 



g s @m 



X- 



t sin 



xf!(t) 



2 



cosh 



g s $rn 



x. 



m 



t sin 



/ m<pt 
V 2 



— sinh 



g s $m 



x. 



t cos 



(61) 



(62) 



Here we have taken the mode function at the initial time as a free field. We can show that 
these functions satisfy the exact scaling property (l59l) . See the details in Ref . [T7I) . Moreover, 
substituting Eqs. flBT]) and (1621) into Eq. fl60|) . we obtain the occupation number at t — nT 
with exact scaling property as 



f x ( n T, K) = sinh 2 



9s&rn x 



nT 



(63) 



This result is consistent with Eq. (157)) . Therefore, the analytical result for the occupation 
number obtained by the time averaging method can satisfy the exact relation in its evolution 
coming from the periodicity of the <fi oscillation. 
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4.2. Growing mode of fermion 

Next, we turn to consider the occupation number of ip with k = at later time. The 
equation of motion for the mode function is now given by 

d 2 P(r) 



dr 2 



+ 



1 + 2 J I - 0* cos(2r) - 2i sin(2r) + g 2 cos 2 (2r) 



P(r) 







(64) 



where q^ = 2g F ^>/m^. As before, we shall use the method of variation of parameters and 
write 



P(t) = Z/i (t) Ui(r) + 2/ 2 (r) « 2 (t) 



(65) 



where y\ = cost and w 2 = sinr are the solutions of (164)) for the case = 0. Together with 
the condition fl53|) we obtain the equations for u\ and u 2 as 



d I U\ 
dr I W2 



2g^, 



/3j cos(2r) - 2i qi , sin(2r) + gJcos 2 (2 



SUIT COS T 



— cos r 



sin 2 r 



SUIT COS T 



Mi 

U 2 



(66) 



By integrating over the period of the <fi oscillation, the averaged U\ and w 2 satisfy 



d 1 






aW \ 








(67) 



From the initial conditions (132 j) at the leading order, namely, P(0) = (1 + (1 — Z? 2 ) 1 ^ 2 ) 1 ^ 2 
and dP(0)/dr = —iP(0), we obtain the solution 



P(0) 



+ i 



1-^-1 



e-sin f^M 
V 2 



(6? 



It is then found from (1341) that the occupation number for the growing mode is given by 



fy(t, K) = sin 



■2 I QtpPi> T 



. 2 ( g F $^ 
sin I — - 1 



(69) 



Note again that this reproduces the initial behaviour in Eq. ( I4"8"j) for t <^ t CT = 2/ (g F ^>(3^), 
as in the scalar production. In Fig. [SJ z CT is around 300 when g F = 10~ 8 and (3^ = 0.5. On 
the other hand, for t 3> t CI , the occupation number oscillates around = 1/2 and does not 
exceed one, which should be contrast to the scalar production. This behaviour reflects the 
Pauli blocking effect of the produced fermion ip. As a result, the number density of ip stops 
to grow at t ~ t cr even if the coupling g F is extremely small. 
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As in the case of scalar production, it can be confirmed that the solution in Eq. (|68|) 
satisfies the exact relation between the mode functions at different times. For all the mode 
P(t, k) can be written as 



P(t,k) = P(0,k) [P (1) (t,k) -too k (0)P {2) {t,k)] . 



(70) 



Here P^' 2 \t,k) correspond to linear independent solutions of the equation of motion with 
the initial conditions pW(0,Jfe) = pW(0,k) = 1 and pV>(0,k) = P {2 \<d,k) = 0. These 
solutions satisfy the same scaling property as Eq. fl59|) by replacing xi -P^ 1 ' 2 ^ and the 

occupation number of ip at t = nT is obtained by same manner as 

k 2 sin 2 (nd) 



f^ipT, k) 



co k (ry S in 2 (d) 

sin(n d) 



[lmPW{T,k)] 
U({n-l)T,k) 



(71) 



sin ((n — 1) d) 

where d = cos^ 1 (ReP^ 1 -* (T)) . This equation has already been presented in Ref . IT2]) . however, 
the estimation of P^(T, k) is essential to obtain f^(nT, k) by using Eq. (ITT)) . 

Now the time averaging method gives us the analytical expression of P(t, fc*) for the 
growing mode. Therefore, we can estimate f^(nT, fc*) analytically by taking P^'(t, fc*) and 
P (2) (t,K) as 



P^ ^ (t, k* 



cos 



cos 



i— sin 

Pip 



sin 



P (2 \t,K) 



2 

m<t> Pip 



9f & Pip , 
sm -t 



(72) 



cos 



sm 



H sin 



cos 



9f$< 



9F$PiP_ 



tl+ ^ Sin 



9F$Pip 



(73) 



Substituting Eq. (172|) into Eq. ( ITT]) , we obtain the occupation number at t = nT with exact 
scaling property as 



U(nT, k*) = sin 2 



2 / g F $/3ip 



nT 



(74) 

obtained by the time 



Similar to the scalar production this result is consistent with 
averaging method. 

Thus, we conclude that the results by the time averaging method can not only abstract 
the characteristic evolution due to the non-perturbative correction but also satisfy the exact 
scaling property in terms of oscillation period of </>. 
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§5. Conclusion 



We have investigated the particle production from the coherent oscillation by using the 
method based on the Bogolyubov transformation. For the case when the coupling constants 
of the oscillating field are very small, we have obtained the leading contributions to the 
distribution functions and the number densities of the produced particles. 

When the amplitude of the oscillation is small (<P <C (</>)), the leading contributions to the 
yields are found to be 0(g 2 s F ). We have presented the exact expressions for the distribution 
functions of the produced x an d ip at the 0(g 2 SF ) order. It has been shown that there exists 
the growing mode with k — k* if vn x ^ < m^/2, and its occupation number increases at 
the rate t 2 after sufficient numbers of the oscillation. The distribution function has a peak 
at k ~ k* and the width of the peak decreases at the rate 1/t. As a result, the number 
density of produced particles is proportional to t. The expression for the number density is 
found to be consistent with the one obtained by assuming that the coherent oscillation is a 
correction of non-relativistic scalar particles and the decay process is a main source of the 
particle production. 

We have found that the above perturbative results fail to describe the exact ones for 
sufficient late times since the non-perturbative correction becomes significant even when the 
coupling constants are extremely small. Indeed, the occupation number of x f° r the growing 
mode increases exponentially while that of ip oscillates around 1/2. These distinctive features 
represent the statistical properties of the produced particles, i.e., the effects of the Bose 
condensation for the scalar production or the Pauli blocking for the fermion production. Due 
to these non-perturbative effects, the explosive production of x happens while the production 
of ip becomes insignificant for late times. To handle with these non-perturbative effects we 
have used the time averaging method, and have successfully described the evolution of the 
occupation number for the growing mode. This method works well because the typical time 
scale of the evolution is much longer than the rapid <ft oscillation. Furthermore, we have 
shown that the results obtained by the time averaging method satisfies the exact scaling 
properties in Ref JT6j) . which also gives the justification of the use of the time averaging 
method. 

Throughout this analysis, we have neglected the back-reaction effect of the produced 
particles in the estimation of the yields. When the occupation number of these particles is 
close to unity, such an effect should be taken into account. In addition to this, the inclusion 
of the expansion of the universe is also necessary to reveal the reheating/preheating processes 
in the inflationary universe. These issue will be discussed in elsewhere.^ 1 
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Appendix A 

Derivations of Eqs. /[2b]) and / f50j) 



We show here the derivations of the number densities n x and at the leading order 
given in Eqs. fTSBI) and floUl) . Let us first consider Eq. fl2"B"|) . It is found from Eqs. f[T$j) and 
(I2T1) that the leading 0(g 2 s ) contribution to n x is given by 



nJt) 



g 2 s <P 2 m\ m\ 



KJt) , 



where 



K x {t) = I dti sin(m0ti) / dt 2 smfa^) J x (At) , 
Jo Jo 

f°° k 2 
J x {At) = dk — cos(2w x At) , 

JO U Y 



with At — t-i — 1\. The integration in Eq. (IA-3[) can be done as 

71 



JJAt) 



4m. 



14 \m'i J/' 2 - I ni x \Al\ ,/ \ ( --;l,-;-m x At 2 



(A-l) 

(A-2) 
(A-3) 

(A-4) 



where 1-^2(0; b, c; x) is the generalized hypergeometric function. We then expand J x in terms 



of m x as 



7T 



J x(^) = — 



- + (Z\tm x ) 2 -|Z\t|m x ^ 



(_l)n+i 



n=0 



(4n 2 -l)(n!) : 



-(Atm x ) 2n 



(A-5) 



Now, we can perform the integrations of ti and £2 m (IA 2[) as 

Jf x (f) = * t [l - 2r 2 x - 2r x - 4r x - 10r x + 0(rf)] + 



7T 

2 



■t& + ... 



(A-6) 



where r x = m x /m ( j ) . Notice that we we have only listed the terms proportional to t and 
dropped off the oscillation terms. Finally, we obtain the number density up to 0(g 2 s ) as 
Eq. 026]) 



n x {t) 



ql <P 2 m 2 m\ ql $ 2 m 2 B y 

x *-K x (t)~- x ' 



8tt 2 



8tt 



(A-7) 
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by neglecting the oscillation terms. 

Next, we turn to consider Eq. fl50|) . As in the case of the scalar production, the leading 
O(gp) contribution to the number density of ip is given by 



ql $ 2 ml 



(A- 



where a prefactor 4 counts the internal degrees of freedom of ip and 

Kjf,(t) = / dtx sin(m^i) / dt 2 sm(m^t 2 ) J^At) , 
Jo Jo 

f°° k 4 
J^(At) = / dk — cos(2uj^At) . 

Jo ^4, 



(A-9) 
(A-10) 



Notice that, comparing with Eq. flA-4[) . the integrand of has an extra factor k 2 . 

When m^j = 0, we first integrate ti and t 2 in and then the k integration in gives 



(All] 



On the other hand, when ^ 0, we first estimate J^(At) as 



M At ) 



71171^ 

4 

nm^ 
4 



3 + 4m^Z\t 2 - 6m^|^|iF 2 



3 + Arn^ At 2 ~ 6r M^*l 5^ 



— ; -.2; -m 2 ,,At 2 
2 2 * 



n+l 



n=0 



(4n 2 -l)(n + l)(n!) 2 



(Atm^) 



2n 



(A-12) 



After the integrations over ti and t 2 we find apart from the oscillation terms 



= ^ l-K + 3rJ + 4rJ + 3rJ + (9(4°)] + . . . , 

where = m^jm^. Combining the above two cases, we obtain 

K4t) = ^ [1 - 6rJ + 6rJ + 4rJ + 6rJ + 0{rf)] +... 
- \t$l + • • • , 



(A-13) 



(A- 14) 



which gives Eq. (150 



?v(t) 



g 2 F P\<Z> 2 m\ 



(A-15) 



by neglecting the oscillation terms. 
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